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1. Introduction 



The purpose of this paper is to describe two methods 
for constructing smooth bivariate functions which take on 
given values at scattered points in the plane. Given the 
data (x^, y^,f^), i = 1,...,1J, both of these schemes 
define a smooth bivariate interpolant S with the property 
that S(x^,y^) = f^, i = 1,...,N. 

In the past several years a number of methods have been 
proposed (e.g. , [1], [3], [6], [7], [8], and [11]), and two 

survey papers, [2] and [10] , have dealt extensively with this 
topic. The problem of constructing smooth approximations 
based upon scattered data is encountered frequently in many 
areas of scientific applications. Common examples are: 
meteorological information such as rainfall and solar 
insolation, geographical information such as elevations, 
geological information such as depths of underground 
formations, and engineering data such as stress values 
obtained by finite element analysis. A somewhat less obvious 
example is given in [10] where it is described how human 
heart potential is measured at irregularly spaced points as 
an aid in diagnosing abnormal heart conditions. 

Since a number of methods are available for this 
important problem, a project was undertaken to test and 
compare as many methods as possible [4] . While the project 
included both global methods (meaning that S(x,y) is depend- 
ent on all data points regardless of their distance from 



(x,y)) and local methods (S(x,y) does not depend on data 
points more than a certain distance from (x,y)), for large 
sets of data it is necessary to use local methods. During 
the course of developing and testing variations of previously 
published schemes we have discovered two which appear to be 
preferable as general purpose methods. While neither of 
these methods have appeared in the literature, they are both 
modifications of previously described techniques. In Section 
2, we describe the basic method from which both of our methods 
derive. In Section 3, we describe the details of the two 
modifications and in Section 4 we show the results of applying 
these methods to certain test problems. In Section 5, we 
discuss some generalizations and ways to fine tune the 
schemes to suit the particular needs of a user. 

2. Inverse Distance Weighted Least Squares Interpolation 
It is convenient to consider the data as coming from 
an underlying function f and to view the interpolation 
process as an operator applied to this function. That is. 



p^(x,y) denote a function with certain properties of a 
distance function. In particular, we assume that p^(x^,yjj^) = 
0 and that l/p^(x^,y^) is a nonnegative decreasing function 
as (x,y) "gets further" from (x^,y^). For example, the 



is one possibility for p^. We denote by ^ , j = l,...,m 
a set of basis functions to be used for a least squares 
approximation . 



fi = f (x^,y^) , i = 1 



N and S(x,y) = P[f](x,y). Let 



usual Euclidean distance, d^(x,y) 




2 



3 

McLain [7] was the first to consider a family of inverse 
distance weighted least squares approximations. The general 
form of these interpolants is 



m 



(2.1) P[f] (x,y) = I a . (x,y) ()) . (x,y) 

j = l ^ ^ 



where aj(x,y), j = l,...,m represents the solution of 



(2.2) Min 



N 

I 



af, . . . ,am i-1 



a,(b,(x.,y.) +...+ a (b (x.,y.) - f. 
1^1 1 -*^1 m^m 1 



(x,y) 



for given (x,y) . McLain gave results of a number of tests 
for several choices of basis functions, (t)^, and distance 

functions, p.. The <}) . consisted of the low order monomials, 

k 'C- 2 ctd 

X y , and was taken as d^, d^ , or d^e i . The higher 

power of d^ and the exponential are motivated by the desire 

to place less importance on distant data than would be 

accomplished by d^ alone. 

The function (2.1) is not defined at any data point, 

(Xi,yi), by the above definition. Because the weight for 

2 

the best least squares approximation 1/p^ (x,y) — ► «> as (x,y) — »■ 
(Xi,yi) it is clear that P[f](x,y) — f ^ as (x,y) — ► (x^,y^), 
else the sum of the squared errors would be unbounded. Thus 
we obtain a continuous approximation if we define P[f] (x^,y^) = 
f^, i = 1,...,N. McLain asserts that the interpolants are 
infinitely differentiable, and upon the assumption of non- 
singularity of the coefficient matrix for the normal equations 



obtained from (2.2), this readily follows- 

McLain singled out the case of the bivariate quadratic, 

P[f] (x,y) = a^^ + a 2 X + a^Y ^s^Y ^6^^ ' with 

ctd ^ 

weight function = d^e i as working well for a variety 
of problems in the sense that small deviations from certain 
test functions were observed. Additional experimentation [3] 
has confirmed this, but even so this method has two serious 
drawbacks: (1) The computational effort required is large 

since each evaluation requires the solution of a least squares 
problem, and (2) the method is global, that is, the interpolant 
depends on all data points regardless of how far away these 
points are from the point at which the interpolant is being 
evaluated. 

Both of the methods we propose can be viewed as modifica- 
tions of the above inverse distance weighted least squares 
interpolant in that they consist of replacing aj(x,y) with 
an approximation A[aj] (x,y) which is computationally more 
tractable than a^ (x,y) itself. We note that as long as 

A[aj ] (x^,y^) = (x^,y^) , i = 1, . . . ,N, 

the operator 

m _ 

Q[f] (x,y) = I A[a .] (x,y) (|> . (x,y) 
j = l ^ ^ 



will maintain the interpolatory properties of P[f]. The type 



of approximation we use for a ^ , j = 1 



,m is of the form 



_ N _ 

A[a.] (x,y) = I a. (x. ,y.) W. (x,y) 

J i=l J 1 1 1 



where 



(2.3) 



W.{x.,y.) = 6.., i, j = 
i j j ij 



Therefore, we may rewrite the above expression for Q as 



N 

(2.4) Q[f] (x,y) = I W. (x,y)Q. (x,y) 

i=l ^ ^ 



where Qj^{x,y) 



m 



I 

j=l 



(x^,y^) (j)j (x,y) 



is the inverse distance 



weighted least squares fit at the point (x^,y^) . 

We refer to the functions Q^, i = as the nodal 

functions since they are associated with the nodes (x^,y^) , 
i = 1,...,N, respectively. We note that Qj^(x^,y^) = f^, 
i = 1,...,N, and that Qj^(x,y) is a local approximation (near 
(Xi,yi)) to f(x,y), and as such may be expected to mimic 
the shape of f(x,y) provided distant points do not influence 
Qj^(X/Y) too much. 

In order for the interpolant Q[f] to maintain the local 
shape characteristics of the nodal functions we will require 
certain properties for the W^ in addition to (2.3) . Specif- 
ically, to preserve 
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9Q[f] 

3x 







3Qi 



and 



8Q[f] 

3y 



(x.,y.) 



3Qi 

3^(x..y.), i . 1 



,N 



we require that 



(2.5) 



8W. 8W. 

ir (x.,y.) = ^(x.,y.) = 0, 1, j = 1 , . . . ,N . 

8x j j 8y j j ' ' J 



For our two methods, we also propose the use of bivariate 
quadratics for Qj^(x,y), i = Thus, if f itself is 

a quadratic function, the function will be identical to 
f, i.e., Q^(x,y) = f(x,y),i = 1,...,N. 

Therefore 



N N 

Q[f] (x,y) = I W. (x,y)Q. (x,y) = f(x,y) I W. (x,y) 
i=l ^ ^ i=l ^ 



whenever f is a quadratic function. In order to obtain 
quadratic precision of the modified interpolant Q[f], we add 
the additional requirement that 



N 

(2.6) I W. (x,y) E 1. 

i=l ^ 



The choice of the distance function to be used in 
(2.2) when calculating the nodal functions was made on the basis 
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of extensive numerical testing [4]. Franke and Little [2, 
p. 112] proposed 



(2.7) 



1 . - 'il)t 

Pi 




R. 



1 



R. 



1 



- d. > 0 

1 

- d^ < 0 , 



and for appropriate choice of the values R^ , this works well. 

We recall that is the solution of the inverse distance weighted 
least squares problem at (x,y) = • Discussion of the 

is simplified (as is the actual computation) by assuming 
that 



Q|^(x,y) = 






a, . (x-x, ) 
k4 k 



ak 5 (x-x^) (y-Y^) + ^k6<y-yk>^ 

We then seek the solution of 



Min 

^k2' • * • '^k6 



N 

I 

i=l 



^k + ^k2(^i-^k^ ^ke^^i'^k^ - ^i 



Pi(Xk,Yk) 



If — is given by (2.7), then whenever d. (x, ,y. ) > R. , 

P • Q. K. K. J. 

the point (x^,y^,f^) has no influence since the corresponding 
term in the sum is zero. Thus depends only on "nearby" 
points and is therefore a local approximation to f. If we now 
use weight functions W^^ which are nonzero only in some neigh- 
borhood of (^j.'Yj.) we will obtain a local interpolant. 

The proper choice of the "radii of influence", R^, 
is critical to the success of the method, and we will discuss 



8 



this in the context of our first method in the next section, 
as well as in Section 5. 



3. Two Methods for Interpolation To Scattered Data 

Specifying the weight functions W^(x,y), i = 1,...,N to 
be used in (2.4) will define an interpolant. We have found 
that both of the choices we propose have very comparable 
fitting capabilities, but we feel that there are situations 
in which one or the other may be preferable, and so both are 
presented . 



3 . 1 Method I . 

This scheme is based upon a special case of the inverse 
distance weighted least squares interpolant given by (2.1) - 
(2.2). If m = 1 and 4)j^(x) = 1, then (2.2) can be explicitly 
solved to yield 

N f . 

i=l p? (x,y) 

' -fj— T 

i=l (x,y) 



and so 



(3.1) p[f] 



H 



I 

i=l 



N 



y 

i=l 
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This method was first proposed by Shepard [11] . Without 
modification it does not have very good fitting properties. 

Gordon and Wixom [5] have analyzed this method and have proposed 
some interesting modifications. They also discuss some application 
areas that are well suited for this type of interpolant. The 
fact that Shepard's method is a special case of inverse distance 
weighted least squares interpolants has been pointed out by 
several authors, e.g. [10]. It can easily be verified that as 
long as p^(x^,y^) = 0, the functions 

(3.2) W. = Pj 

k=l pj 

will satisfy the conditions of (2.3), (2.5), and (2.6). This 

last statement assumes that the distance functions p^, i = 1,...,N 
are sufficiently smooth so that the derivatives (2.5) exist. 

This is certainly the case for our choice of p^, i = 1,...,N 
given by (2 . 7) . 

We complete the description of our first method with a 

discussion of the selection of the involved in the 

definition of p^. While the use of variable radii (i.e., 

Rj) adds to the flexibility of the methods, we have found that 

for a general purpose interpolant, the selection of these values 

can quite often be a nuisance which can be avoided by the use of 

a uniform value (R. = R for all i) . 

1 

At this point we emphasize that the distance functions 
p^ enters in two places: (1) definition of the nodal 



10 



functions, Qj^, and (2) definition of the weight functions, 

It is not necessary to use the same radius of influence for 

both instances, and experience has shown it is desirable to use 

different values, say R = R in defining the nodal functions, 

and R = R^ in defining the weight functions. Since R^ denotes 

the radius of influence of the data points on the nodal functions, 

while R denotes the radius of influence of the nodal functions 
w 

on the interpolant Q[f] (x,y) , it is clear we should take R^ < 

Rg. In order to aid the naive user in making reasonable choices 
of R and R we have found it convenient to specify values of 
Ng and and to compute the radii of the influence regions 
according to the relationships: 






(3.3) 



w 2 V N 



where D = 



max d^ (x . ,y . ) . 

i/ j ^ 



These choices of R and R eliminate the effects of scaling 

q w ^ 

the data. The values of N and N can be thought of as 

q w ^ 

representing the number of data points anticipated to lie in 
circles of radii R^ and R^, respectively. For somewhat uniformly 
distributed data, we have found that a value of N =18 works 

q 

quite well. For data that have some regions which are 
relatively sparsely populated and other regions where the 
data are comparatively dense, or for small sets (N < 25), it 
may be necessary to increase the values, since the interpolant 
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is defined only on the union of disks of radius R centered 

at the data points i = 1,...,N. We have also found 

that the use of the relationship N /N 2 is useful. To 

q w 

avoid problems when fewer than five additional data points 
fall within a distance of some given (x^,y^) , we have 
incorporated an automatic fallback to linear nodal functions 
in this case. In general, we use the singular value decomposi- 
tion to compute the coefficients of the nodal functions, which 
avoids possible nonuniqueness problems. It is not usual for 
either situation to occur, however. 

We now summarize the description of our first method. 

i) Select N and N in order to define 
q w 




(R -d.)^ 
R d. 

q 1 





(R -d.)^ 
w 1 + 

R d. 
w 1 



R 



w 




Default values of and are 18 and 9, respectively, 
ii) For k = 1,...,N solve the least squares problem: 



N 



min 



•) , 



i=l 









-I 



■"®k4 ^+®k5 <’‘l-^k> 'i’i-^k’ +®k6 ‘yi-^k' 
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to yield ^ , j = 
iii) Define 



2 , 



6 . 



Q^(x,y) = + ^k2<=<-=<k> ®k3<y-yk> 

+ ir,^^(x-x^)2 + S;^5(x-x^) (y-y^) + 

and compute 



D[f] (x,y) = 



N Qj^(x,y) 

kliP^T^ 



N 

I 

k=l 



111 Pv (x/y) 



3.2 Method II. 

Our second method makes use of a triangulation of the 
data = (x^,y^) i = 1, . . . ,N in order to define the weight 
functions W^, i = 1,...,N. We use an algorithm which triangulates 
the convex hull based upon the min-max angle criterion as 
described by Lawson [6]. A FORTRAN program which implements 
this algorithm is available as part of [1]. Alternatively, if 
a triangulation is in existence for other purposes, it can be 
used. 

Each will be a globally defined function with support 

S. = U T..,, where T.,, denotes the triangle with 

jkl e 

vertices V., V, , and V, and M. = {jkl: T 

3 K 11 

with vertex V^}. 



is a triangle 
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Figure 1: 



We first define and its first order partial derivatives on 

E = U e, . where e, . represents the edge with vertices 
kjeN 

e 



V, and V. and N = {kj: V, V . is an edge of the triangulation}. 

K J © JC ^ 

Following this, we incorporate a blending method for triangles 
to extend the definition to the interior of each triangle of 
the triangulation. Let e^^ be an edge contained in with 
as an endpoint. As a univariate function along this edge, 
must satisfy four conditions imposed by (2.3) and (2.5). 
Namely 



W^(V^) = 1, W^(Vj) = 0, 



aw. aw. 

(x.-x. )ir— ^^(V. ) + (y .-y . )-r— ^(V. ) 

3 1 ax 1 -^3 ■'i ay 1 



0 



and 




aw. 



aw . 



0. 
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These conditions can be satisfied by a cubic polynomial and 
so we define 

W^((l-t)V^ + tVj) = (1-t)^ (2t+l), 0 < t < 1. 

On all other edges we define to be zero. In order to 

maintain continuity of the first order derivatives across 
edges, it is convenient to specify the derivatives normal to 
an edge as a linear function along the edge. That is 

8W. 9W. 

3W. 3W. 

= (l“t) [ (y . -y . ) — ^(V.) - (x.-x.)-r — ^(V.)] 

■^j -^1 9x 1 j 1 9y 1 

9W. 9W. 

+ t[(y.-y.)- 5 — ^(V.) - (x . -X . ) V (V . ) ] ♦ 

■'x 9x ' 3 i' 9y j' ^ 



In light of (2.5), this means that the normal derivatives will 
all be identically zero. This completes the description of 
the edge information for W^. 

In order to extend the definition of W. to the interior 

1 

of each triangle, we use an interpolation method [9] which 
will assume predescribed position and slope on the entire 
boundary of a triangle domain. After substituting the edge 
information into this triangular blending method, we find 
that for (x,y) e T. c s., the weight functions have the 

1 j K 1 



VJ. (x,y) = b. (3-2b. ) 
1 1 1 



form 
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(3.4) 



+ 3 



b?b .b, 
13k 



b . b . +b . b, +b . b, 
13 1 k 3 k 



b . 
3 



^k 







21 




where b., b., b, are the barycentric coordinates of (x,y) 

1 3 K 

with respect to the triangle and ll^j^ll n = i, j or k 

represents the length of the edge opposite V^, n = i, j or k. 
The barycentric (area) coordinates are given by the equations 





X = b . X . 


+ b .X . + 


b, X, 




1 1 


3 3 


k k 


(3.5) 


Y = b . y . 
^ i-' 1 


+ b . y . + 
3-^3 






1 = b. + 
1 


b . + b, 

3 k 


• 



We can now note that on an arbitrary triangle the only 

weights which are nonzero are W^, and W^. Therefore, the 
final interpolant is given by 



(3.6) G[f] (x,y) = W^(x,y)Q^(x,y) + (x , y ) (x , y) 

+ (x,y) (x,y) , (x,y) e 

We also note that W. + W. + W, = 1, so that (2.6) is satisfied. 

1 3 K 

We now summarize this method. 

1) Define the nodal functions Qj^, k = 1,...,N 
as in Method I, 

ii) Form a triangulation of the points = (x^,y^) , 

1 1,...,W, 

Given (x,y) determine the vertices V., V. and V, 

13 K 

of the triangle that contains this point and 
compute G[f](x,y) according to (3.6) using (3.4) and (3.5). 



iii) 
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We note that this method is very similar to that proposed 
by McLain [8]. Our weight functions W^, i = 1,...,N arise in 
a natural manner, however, and use of the distance function 
given in (2.7) is an improvement since it is generally the case 
that shape characteristics of the nodal functions in relation 
to f are adversely influenced when global approximations 
are used. 



4 . Examples 

In order to illustrate the performance of the two methods, 
we include some examples. These examples are only a few of 
many upon which our conclusions are based, but are representative 
of the methods' approximation properties. The first group of 
examples utilizes ordinates from the function 



f (x,y) 



. 75EXP 1^- - I - j + .75EXP 

.2EXP f-(9x-4)^ - (9y-7)^j+ . 5EXP 



(9x+l)^ (9y+l) 1 

49 10 J 



(9x-7)^+(9y-3) 

4 




A perspective plot of this surface, viewed from the first 
quadrant at an angle of 30° from the x-axis is given in figures 
3a and 4a. The function was chosen so as to present a variety 
of behavior in a single surface. The maximum value of the 
function is approximately 1.22 near the point (.22, .22), while 

the minimum value is approximately .004 near the point (.47, .78). 

Three sets of data (x^,y^), i = 1,...,N were used. 

Set 1: This set consists of 100 points generated by a 

pseudorandom number generator, one point in each subsquare 
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(a) Data Set 1 
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(b) Data Set 2 



(c) Data Set 3 



Figure 2 



of side 1/9 centered at (i/9,j/9), i, j = 0,1,.. .,9. These 
points are shown in Figure 2a. 

Set 2: This set of 33 points was selected manually to have 

regions of varying density. These points are shown in 
Figure 2b. 

Set 3: This set of 25 points was selected manually to yield 

a somewhat uniform coverage of the unit square, and are similar 
to a set appearing in [8] . These points are shown in Figure 2c 

The interpolants were evaluated on a uniform grid of 
33 X 33 points in the unit square. The resulting surfaces 
are shown in Figures 3 and 4. Table 1 contains the maximum, 
mean, and RMS deviations at these points. For Method II a 
few of the display points lie outside the convex hull. For 
plotting purposes, these were set to zero and were omitted 
in the calculations leading to the errors of Table 1. 



Case 


Max 


Mean 


RMS 


1— 1 

• 

H 


. 0573 


.0079 


.0128 


CN 

• 

H 


.1844 


.0340 


.0478 


1.3 


.1584 


.0353 


.0486 


II. 1 


.0481 


.0072 


.0113 


II. 2 


.1501 


.0326 


.0455 


II. 3 


.1535 


.0349 


.0475 



Table 1. Errors 
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Figures 3b and 4b show that both methods reproduce the 
surface quite well for data set 1. The main defect appears 
to be near the higher peak, particularly for Method I. As 
can be seen from the disposition of the points in Figure 2a, 
this is in a region where a relatively large gap between points 
occurs. A similar occurrence accounts for the depressed area 
behind the dip in both figures 3b and 4b. 

Figures 3c and 4c are quite similar and both have noticeable 
defects when compared to the test surface. In particular, the 
dip is completely missed because the data points fail to define 
it. The appearance of a crease in Figure 4c near the right 
rear edge is due to the occurrence of a long thin triangle 
along that edge which causes blending of two nodal functions 
from points relatively far apart in the definition of the 
interpolant for that triangle. When these two nodal functions 
differ significantly, as they do here, rapid variations can 
occur across the narrow part of the triangle. 

Figures 3d and 4d appear to be less alike than they actually 
are because values outside the convex hull have been set to 
zero in Figure 4d. The most significant difference between 
the two is near the left rear edge, where Figure 3d shows 
the surface (apparently) dipping down rather rapidly, while 
Figure 4d shows a near crease similar to that in Figure 4c. 
Because of a data point near (.47, .78), the dip is partially 
defined in this case, but since there are no other nearby 
points it is extended over a much wider area than in the 



test surface. 
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Figure 3: Method I 
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(a) Test surface 




(b) Interpolant for Data Set 1 




c) Interpolant for Data Set 2 (d) Interpolant for Data Set 3 



Figure 4 : Method II 



Another set of data obtained from Akima [1] , which arises 



in a study of waveform distortion, is given in Table 2. We 
use this set to illustrate the effects of varying the parameters 
N and N . The results are shown in Figure 5. 

Figures 5a and 5c appear to be difficult to choose between. 
Very slight differences can be observed along the front edge. 
Figure 5e is definitely less desirable than 5a or 5c because of 
more undulations near the front edge and a higher peak at the 
right rear. Extensive testing has shown that Method I is 
fairly stable for values of and in the ranges given here. 

Figure 5d appears to be the more desirable surface among 
5b, 5d, and 5f. Figure 5b shows a small defect near the right 
front edge, while the choice between 5d and 5f is less obvious. 
All three show the characteristic defect over the long slim 
triangle at the right rear edge, allowing the surface to dip 
out of sight there. It should be pointed out that there is 
no known "parent" surface here, and in fact that behavior may 
be correct, although that part of the interpolant is not 
pleasing because of the very rapid changes which occur. 



5. Conclusions and Recommendations 

The two schemes discussed here have been found to be 
capable of generating reasonable interpolation functions in 
a variety of cases. A number of comments are appropriate, 
however . 

First, we have restricted ourselves here to the 
discussion of local interpolants . Local methods are necessary 
for very large sets of data, but in exploring their properties 



i 


X . 
1 




z . 
1 


i 


X . 
1 




z . 

1 


1 


11.16 


1.24 


22.15 


26 


3.22 . 


16.78 


39.93 


2 


24 . 20 


16.23 


2 .83 


27 


0.00 


0.00 


58.20 


3 


12.85 


3 . 06 


22 . 11 


28 


9.66 


20.00 


4.73 


4 


19.85 


10.72 


7 .97 


29 


2.56 


3 . 02 


50.55 


5 


10.35 


4 . 11 


22.33 


30 


5.22 


14 . 66 


40.36 


6 


24.67 


2.40 


10.25 


31 


11.77 


10.47 


13.62 


7 


19.72 


1.39 


16.83 


32 


17 . 25 


19.57 


6.43 


8 


15 . 91 


7.74 


15.30 


33 


15.10 


17 . 19 


12.57 


9 


0 . 00 


20 . 00 


34 . 60 


34 


25 . 00 


3.87 


8.74 


10 


20.87 


20.00 


5.74 


35 


12.13 


10.79 


13.71 


11 


6. 71 


6 . 27 


30.97 


36 


25.00 


0.00 


12.00 


12 


3.45 


12.78 


41.24 


37 


22 . 33 


6.21 


10.25 


13 


19.99 


4 . 62 


14.72 


38 


11.52 


8.53 


15.74 


14 


14.26 


17 . 87 


10.74 


39 


14 . 59 


8.71 


14 . 81 


15 


10.28 


15.16 


21.59 


40 


15.20 


0.00 


21.60 


16 


4 .51 


20 . 00 


15 . 61 


41 


7.54 


10.69 


19.31 


17 


17.43 


3.46 


18 . 60 


42 


5.23 


10.72 


26.50 


18 


22 . 80 


12.39 


5.47 


43 


17.32 


13.78 


12.11 


19 


0.00 


4 . 48 


61.77 


44 


2.14 


15.03 


53.10 


20 


7.50 


1.98 


29.87 


45 - 


0.51 


8.37 


49.43 


21 


16.70 


19 r 65 


6.31 


46 


22.69 


19.63 


3.25 


22 


6. 08 


4 . 58 


35.74 


47 


25.00 


20 . 00 


0.60 


23 


1 . 99 


5.60 


51.81 


48 


5.47 


17.13 


28.63 


24 


25 . 00 


11.87 


4 .40 


49 


21.67 


14.36 


5. 52 


25 


14 . 90 


3 . 12 


21.70 


50 


3.31 


0.13 


44.08 



Table 2 . 




(a) Method I: N =12, N =6 
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(b) Method II; N =12 
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(c) Method I; N =18, N =9 
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(d) Method II: N =18 



Figure 5. 
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Method I: N =24, N =12 
q w 




(f) Method II 



N =24 
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Figure 5. (continued) 



one need not consider large sets because they are local. 

We have only considered interpolation methods here, although 
we recognize that often it may be more appropriate to use 
approximation methods which smooth the data in some sense. 

While we have not investigated this possibility, it is clear 
that replacement of the nodal functions Qj^(x,y) with local 
smoothing functions rather than ones which take on the value 
fj^ will lead to a smoothing approximation which is local. 

The use of the same radius of influence for each point, 
and the calculation of that radius from the diameter of the 
point set is strictly a matter of convenience for the user. 

One could argue that this device makes the methods global 
since addition of a point will change the radius of influence 
and hence the entire interpolant. For this reason we made the 
computation of the radius of influence an option (although 
it is the default option) . The use of different radii of 
influence could be a necessary and desirable feature when the 
density of points varies drastically over the point set. Our 
experience indicates that one should probably choose radii so 
that the disks associated with a point contain approximately 18 
points for defining the nodal functions (Methods I and II) and 
about 9 points for defining the weight functions (Method I) . 

Regarding the choice between Method I and Method II, 
we make the following comments: Method I has the advantage of 

simplicity. While Method II requires a triangulation (and the 
machinery for generating it if it is not already available) , 



and thus more auxiliary storage, it is considerably faster 
since each evaluation involves only three nodal functions, 
while Method I typically involves about 9 (N^ = 9) nodal 
functions. Method II also has the disadvantages noted in 
the examples when long thin triangles occur. Method II is 
not readily extended to more than two independent variables, 
as is Method I. Nonetheless, for certain applications. Method 
II may be the appropriate choice, particularly if a triangula- 
tion is already in existence. 

In conclusion, we note that all local methods involve 
some ad hoc assiamptions and/or parameters to be specified 
by the user. The methods we propose have endured through 
extensive testing of their fitting properties and appropriate 
values for their parameters. We feel they will perform quite 
adequately in a variety of situations. Nonetheless, we 
recognize that selection of a suitable interpolant is a sub- 
jective matter even in the case of one independent variable, 
and thus the choice of an interpolant ultimately rests with the 
user. 

FORTRAN programs which implement Methods I and II are 
available on request from the authors. 
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